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In  this  work,  we  apply  time-implicit  discontinuous  Galerkin  methods  to  the  problem  of 
thermal  ablation,  where  we  solve  for  the  dynamics  of  flow  of  pyrolysis  gas  in  a  charring 
ablating  media.  We  have  benchmarked  our  results  with  the  published  data.  The  protective 
coating  of  thermal  protection  system  over  the  space  vehicle  is  of  the  order  of  a  few 
centimeters,  and  high  flow  speed  of  gas  and  dominant  source  terms  in  the  problem  put  very 
restrictive  CFL  limit  necessitating  the  use  of  implicit  time  integration  methods.  Good 
convergence  of  these  methods  is  ensured  by  use  of  Newton  method,  where  the  required 
jacobians  are  calculated  numerically  for  accurate  evaluations.  Analytical  approach  to 
evaluate  the  jacobians,  especially  for  thermal  ablation  problems,  is  very  cumbersome,  and 
causes  serious  convergence  issues.  We  have  also  tested  our  development  of  DG  capability  for 
2-D  double  Mach  reflection  problem  for  extending  the  code  to  2-D  thermal  ablation  in 
future. 
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Nomenclature 

maximum  eigenvalue 
speed  of  sound 
specific  heat  of  char 
specific  heat  of  resin 
diffusion  term  for  pyrolysis  gas 
internal  energy  of  char 
internal  energy  of  gas 
internal  energy  of  resin 
void  fraction 
friction  term 

heat  of  formation  of  resin 
Inertia 

thermal  conductivity 
permeability 
eigenvalue 
viscosity  of  gas 
density 


I.  Introduction 


Space  vehicles  upon  entry  into  a  planet’s  atmosphere  typically  experience  high  heat  fluxes  to  their  surface. 
Excessive  heat  flux  results  from  combined  effect  of  convective  and  radiative  components.  Chemical  reactions  of 
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the  dissociated  and  ionized  gas  species  (present  in  the  shock  layer,  which  is  between  the  bow  shock  and  boundary 
layer)  also  add  on  to  the  overall  heating  of  the  surface.  Thus,  role  of  a  protective  coating  on  the  vehicle’s  surface, 
called  Thermal  Protection  System  (TPS)  is  critical  for  the  successful  and  safe  entry  of  the  space  vehicle.  Thermal 
ablation  refers  to  the  mechanism  by  which  ablating  TPS  undergo  degradation  at  high  temperatures,  and  reject  high 
heat-content  material  at  the  surface  itself  thus  shielding  the  vehicle  from  the  effect  of  immense  surface  heating.  In 
addition,  charring  type  ablators  undergo  thermal  decomposition  (of  resin  component)  within  and  release  mixture  of 
gases  called  pyrolysis  gas  into  the  boundary  layer.  These  gases  provide  further  cooling  of  the  medium  both  by 
absorbing  the  heat  from  the  material  and  by  pushing  the  incoming  hot  gases  (from  the  flow  over  the  vehicle)  away 
from  the  surface,  thus  providing  a  reduction  in  overall  convective  heating. 

Problem  of  thermal  ablation  is  crucial  to  the  success  of  TPS  design  studies,  where  the  goal  is  to  accurately 
predict  the  response  of  TPS,  to  given  flow  conditions,  to  arrive  at  an  optimum  mass,  size  and  shape  for  the  TPS. 
Mass  of  TPS  is  directly  linked  to  the  launching  cost  of  the  vehicle,  and  its  size  and  shape  to  the  aerodynamic 
coefficients  of  lift  and  drag.  Thermal  ablation  is  a  very  complex,  highly  non-linear  and  fully  coupled  multi-physics 
problem.  In  addition  to  this,  TPS  for  larger  and  more  advanced  payloads  are  no  longer  a  simple  geometry,  single 
piece  of  material,  but  a  complicated  geometry  and  consists  of  layers  of  stacked  and  multiple  materials.  Most  of  the 
thermal  response  codes  developed  so  far  are  based  on  finite  difference  and  finite  volume  methods  and  are  limited  to 
simple  geometries.  Finite  element  method  has  so  far  been  limited  to  only  conduction  part  of  the  problem,  where  no 
ablation  is  considered.  Also  most  of  the  codes  neglect  the  flow  of  pyrolysis  gas  within  ablating  media  missing  on  its 
effect  on  the  overall  response  of  TPS,  as  shown  in  various  studies.  Capability  of  finite  element  methods  to 
accommodate  complex  geometries  (by  use  of  unstructured  meshes),  ability  to  conveniently  apply  difficult  boundary 
conditions  (typical  of  thermal  ablation  problem  at  the  interface  of  solid  and  fluid  domain)  without  any  loss  of 
accuracy  at  the  boundary,  and  scalability  for  parallel  architectures  (due  to  its  compact  scheme  even  for  higher  order 
schemes)  make  it  particularly  suitable  for  thermal  ablation  problem. 

Discontinuous  Galerkin  (DG)  methods  combine  the  advantages  of  both  finite  element  and  finite  volume 
methods.  They  were  formally  developed  for  hyperbolic  conservation  laws  by  Cockburn  and  Shu  [3-6]  and  later  on 
applied  to  elliptic  or  convection-diffusion  type  of  problems  [7-10].  In  work  by  Persson  and  Peraire,  they  have  been 
applied  to  various  problems  of  viscous  flows,  shocks,  turbulent  flows  and  computational  aeroacoustics  [11-15,  27] 
to  make  them  practical  for  problems  of  practical  interest.  We  extend  the  application  of  discontinuous  Galerkin 
methods  to  the  field  of  thermal  ablation.  As  a  first  step,  we  solve  an  Arc  jet  case  (a  1-D  thermal  ablation  problem) 
given  in  Wakefield  and  Pitts,  1980  [20]  and  later  on  simulated  by  Ahn  and  Park  in  their  two  papers  in  1998  and 
2002  [1,2], 

As  pointed  out  in  [1]  and  [2],  numerical  codes  for  solving  thermal  ablation  problem  assumes  the  flow  of 
pyrolysis  gas  generated  upon  decomposition  of  resin  material  to  be  steady  state.  Such  an  assumption  is  valid  for  high 
heating  rate  environment  and  high  pressure,  wherein  the  gas  leaves  the  system  as  soon  as  it  is  generated.  Thus,  there 
is  no  need  to  solve  for  the  flow  of  pyrolysis  gas  within  the  porous  media.  However  such  an  assumption  is  not  valid 
in  moderate  to  low  heating  environments.  The  governing  equations  of  the  flow  for  such  environment  are  as  given  in 
[2], 
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In  the  above  equations,  Eq.  (1)  defines  the  decomposition  of  resin  material.  Here  pr  is  resin  density  and  R  is 
decomposition  rate  dependent  on  temperature,  T  and  the  resin  density,  pr,  given  in  Eq.  (6).  The  Eq.  (2)  solves  for 
conservation  of  mass  for  gas  density,  here  e  is  void  fraction  and  is  given  by  Eq.  (5).  Here,  £max  is  maximum  void 
fraction  that  happens  in  char  state,  and  pp  is  intrinsic  density  of  resin  and  both  their  values  are  equal  to  0.1788  and 
1763.6  kg/m3.  Opposite  signs  for  R  in  first  two  equations  (1)  and  (2)  indicates  that  resin  decomposes  to  release 
pyrolysis  gas.  Second  term,  D  on  right  hand  side  of  Eq.  (2)  is  the  diffusion  term  that  represents  the  rate  of  change  of 
gas  density  due  to  the  spatial  variation  of  pressure,  and  is  given  by  Eq.  (7).  In  this  expression,  Kg  is  the  permeability 
of  gas,  and  p  is  gas  viscosity.  Eq.  (3)  governs  the  momentum  conservation  of  the  gas.  Flow  of  the  gas  is  pressure 
driven.  Here  u  is  the  gas  velocity  and  spgu  is  the  gas  momentum.  Since  the  gas  is  flowing  through  the  porous  media, 
the  resistance  to  gas  motion  are  given  through  f,  friction  and  I,  Inertia  terms  given  in  equations  (8)  and  (9) 
respectively.  The  Eq.  (4)  solves  for  combined  energy  conservation  of  both  solid  (resin  +  char)  and  gas.  First  term  in 
Eq.  (4)  gives  unsteady  variation  of  total  energy,  E,  of  the  system.  eg,  ec,  and  er  are  internal  energies  of  the  gas,  char 
and  resin  material  respectively.  Second  term  in  Eq.  (4)  is  derivative  of  energy  flux,  similar  to  Euler  equations,  where 
the  flux  is  decided  by  internal  and  kinetic  energy  and  pressure  flow  energy  of  gas.  These  terms  were  neglected  in 
most  thermal  response  codes  since  the  gas  was  supposed  to  leave  the  material  as  soon  as  it  was  generated.  Last  term 
in  energy  equation  is  heat  conduction  within  the  material,  through  which  we  also  apply  net  heat  flux  boundary 
condition  at  the  surface.  The  terms,  ec  and  er,  are  evaluated  using  specific  heats  of  both  materials  (char  and  resin) 
and  are  given  as, 


=  PcCJ 

(10) 

P,C,J  +  PrK 

(11) 

Specific  heat  for  solid  carbon  char,  Cpc  used  in  Eq.  (10)  is  given  as, 


CPC  (r)  =  ^  +  C2^3  +  ^  +  C4Z  +  c5  J-J  (12) 

Here,  z  =  log10T,  Ci  =  -6.5,  c2  =  -143.7,  c3  =  706.7,  c4  =  1835  and  c5  =  -5700.  Specific  heat  for  resin  material,  Cpr 
is  assumed  to  stay  at  constant  value  of  1 174  W/m-K.  This  value  wasn’t  provided  in  [1,  2],  We  consider  two  cases  for 
above  given  governing  equations,  as  done  in  [1]  and  [2],  The  differences  between  both  cases  are  highlighted  in 
Table  1.  Papers,  [1]  and  [2]  use  two  different  models  for  friction  term,  f,  which  have  different  magnitudes.  It  was 
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found  that  the  friction  term  in  [1]  is  two  orders  of  magnitude  higher  than  friction  term  in  [2],  Both  use  two  different 
values  of  thermal  conductivity  and  there’s  no  Inertia  term  considered  in  [1],  Also,  Diffusion  term,  D,  is  missing  in 
the  equations  given  in  [1].  All  these  major  differences  in  two  papers  are  summarized  in  Table  1.  This  table  lists 
friction,  f,  conductivity,  k.  Inertia,  1  and  diffusion  coefficient,  D. 


Table  1:  Differences  in  modeling  terms  in  [  1]  and  [2] 


Difference  # 

1998  paper,  |1] 

2002  paper,  |2] 

1 

f=  7.8  X  107(T/300)°'75(e/0.1788)u 

f  =  2.58  X  105(T/500)'769(e'381)u 

2 

Log10k(T)  =  -6.343X10"yT2  +  3.97X10'4T  -0.193 

Log10k(T)  =  -6.5X10'8T2  +  5X10'4T  -  0.1 

3 

1  =  0 

I  =  (1.222/Ki!0'5)pi;u2 

4 

D  =  0 

D  =  (Ka/|i)(£pB)£02P/dx2 

The  work  in  these  papers  was  a  reasonable  starting  point  for  us  to  test  our  DG  formulation  anchored  in  MIG, 
Multi-scale  Ionized  Gas  flow  Code,  which  is  a  family  of  our  in-house  finite  element  modules,  and  has  been  used  to 
solve  problems  in  plasma  dynamics  like  low  pressure  Hall  thruster,  high  pressure  RF  induced  optical  glow 


discharges,  and  microthrusters,  microchannel  and  nanopore  flows  etc. 
DG  capability  to  use  them  later  on  in  2-D  Thermal  Ablation  problems. 


We  have  also  worked  in  building  a  2-D 


No  flow  outside 
vehicle 


<- 


Q  in  = 
W/cm2 


1400 


Figure  1.  Schematic  of  1-D  thermal  ablation 
problem  (material  is  carbon  phenolic) 


II.  Full  Model  Description 

We  consider  an  Arc-jet  problem  simulated  by 
Wakefield  and  Pitts20  and  Ahn  and  Park1  2.  The  TPS  is  L  (mS'de)  R  (°U'S,de) 

exposed  to  a  heat  flux  of  1400  W/cm2.  The  problem  is  1- 
D  with  domain  thickness  of  1  cm.  As  the  heat  flux 
transfers  into  the  material  through  thermal  conduction, 
the  resin  material  undergoes  decomposition  and 
produces  a  mixture  of  gases,  known  as  pyrolysis  gas.  14 
gas  species  are  considered,  C,  CH,  CH2,  CH3,  CH4,  CO, 

C02,  C2,  C3,  H,  HO,  H20,  O,  and  02.  These  are  taken  to 
be  at  thermal  equilibrium  at  the  temperature  of  the  solid 
material  and  pressure  of  the  gas  inside  the  porous  media. 

Equilibrium  properties  of  the  gas  mixture,  like  molecular 
weight,  internal  energy  and  pressure  are  evaluated  using  an  external  thermo-chemical  solver,  CANTERA21.  The 
properties  are  expressed  as  a  function  of  temperature  and  gas  density.  Further  details  on  thermo-chemical  modeling 
used  in  CANTERA  were  already  discussed  in 
our  earlier  paper,  [24].  We  have  used  two 
curve  fits  for  internal  energy,  eg  (see  Fig.  2). 

First  curve  fit,  (green  curve)  is  a  single 
exponential  curve  for  the  complete  range  of 
temperature  between  200  K  and  3400  K.  This 
curve  under  predicts  internal  energy  for 
temperature  range  from  2400  K  to  3000  K  and 
over  predicts  from  3000  K  to  3400  K.  For 
temperature  range  below  3000  K,  this  under¬ 
prediction  of  internal  energy  is  significant, 
and  results  in  lower  prediction  of  cooling 
provided  by  the  flow  of  gas.  Therefore  for 
more  accurate  prediction  we  pick  second 
curve  as  shown  in  Fig.  2.  This  curve  is  blue 
exponential  curve  from  200  K  to  2800  K  and 
black  exponential  curve  from  2800  K  to  3400 
K.  The  expressions  of  both  internal  energy,  eg 
and  molecular  weight,  MW  are  given  as, 


Figure  2.  Plot  of  internal  energy,  eg,  with  Temperature,  T,  for 
pyrolysis  gas,  comparing  two  curve  fits,  #1  and  #2.  Out  of  both 
curve  fit  -  2  is  more  accurate. 


eg  =exp  [AJ  +  A2)-Aj 


(13) 
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MW  =  B{  exp(jp)  +  i?2  exp(x)  +  53  exp(x)exp(j^)  +  54 


(14) 


Here,  A[  =  1.869e-03,  A2  =  1.198exp+01  and  A3  =  1.156exp+07.  For  Eq.  (14),  x  and  y  depend  on  temperature 
and  gas  density  respectively  as,  x  =  (T/1000)2  and  y  =  LOGi0(pg),  and  constants  are  B!  =  1.7981exp-04,  B2  =  - 
1.333exp-06,  B3  =  2.0159exp-07,  and  B4  =  1.678exp-02. 

Due  to  gas  generation  and  increase  in  temperature,  the  pressure  of  gas  rises  within  the  domain  which  causes  the 
gas  to  leak  out  to  the  atmosphere.  This  gas  as  shown  in  the  results  section  leaks  out  to  the  atmosphere  with  velocities 
of  order  of  100  m/s,  and  within  the  material  domain  they  have  velocities  close  to  20  -  40  m/s,  which  causes 
sufficient  cooling  within  the  material  domain,  while  it  is  ablating.  This  effect  is  neglected  in  other  numerical  works 
in  thermal  ablation. 

III.  Numerical  Scheme  (Discontinuous  Galerkin  Method) 

We  initially  attempted  the  above  described  thermal  ablation  problem  with  Galerkin  based  finite  element  method. 
This  method  develops  spurious  oscillations  in  the  solution  and  is  unable  to  capture  solutions  with  high  gas 
velocities.  Discontinuous  Galerkin  method  allows  for  the  solution  to  be  discontinuous  across  element  boundary  and 
thus  can  accommodate  for  solutions  with  steep  gradients  and  high  gas  velocities.  The  in-built  dissipation  mechanism 
of  DG  methods  results  from  the  use  of  approximate  numerical  fluxes  for  the  inviscid  fluxes.  This  in-built  dissipation 
is  proportional  to  the  jump  in  the  solution  across  solution  variable.  Hence  DG  methods  are  suitable  to  treat  flows 
with  high  speeds.  DG  methods  were  developed  by  Cockburn  and  Shu,  as  a  class  of  explicit  methods,  called  TVD 
Runge-Kutta  Discontinuous  Galerkin  method  (TVD  RKDG),  in  their  series  of  papers  from  1989  -  200 1 3  6.  They 
have  been  applied  to  various  problems  like  Compressible  flow,  Navier  Stokes,  electro-hydrodynamics, 
computational  aeroacoustics,  turbulent  flows,  flapping  wings  etc7'  u"15.  Persson  and  Peraire  have  worked  on  applying 
time -implicit  DG  methods  to  more  practical  and  real  world  problems,  and  some  of  their  work  include  areas  like 
development  of  Compact  Discontinuous  Galerkin  (CDG)  scheme,  preconditioners  for  time-efficient  solution  of 
matrix-systems  obtained  from  application  of  DG  methods,  sub-cell  shock  capturing  by  artificial  dissipation  for 
capturing  flows  with  shock,  and  using  high  order  DG  methods  to  capture  shock  within  one  element  etc  [11  -  15]. 
Their  work  is  very  relevant  to  us  in  the  development  of  implicit  DG  method  for  thermal  ablation  problems. 

A.  1-D  Discontinuous  Galerkin  method 

As  a  first  step,  we  developed  explicit  RKDG  and  DG  methods.  Such  explicit  schemes  are  insufficient  for  thermal 
ablation  problems  as  they  pose  heavy  restrictions  on  CFL  number.  For  example,  for  the  1-D  thermal  ablation, 
maximum  time-step  that  could  be  used  for  the  explicit  DG  methods  was  1.0*  10'08  sec.  Total  run  time  for  thermal 
ablation  problems  is  however  usually  on  the  order  of  few  seconds.  Therefore  explicit  schemes  are  highly  insufficient 
for  solving  thermal  ablation  problems,  even  for  a  1-D  case.  The  situation  worsens  for  higher  dimensions.  The 
severity  of  restriction  on  CFL  numbers  increases  as  the  domain  size  decreases,  and  net  amount  of  gas  flow  velocity 
and  source  terms  increase  in  the  problem.  Implicit  DG  methods  allow  for  higher  time-steps,  e.g.  in  our  current  1-D 
thermal  ablation  problem  we  were  able  to  use  time  step  as  high  as  l.OxlO"03  sec,  which  corresponds  to  a  CFL 
number  of  around  1000.  Since  source  terms  are  very  large,  we  also  get  a  highly  ill  conditioned  matrix,  which  needs 
to  be  solved  using  preconditioners.  A  brief  derivation  of  the  weak  form  of  discontinuous  Galerkin  methods  with 
implicit  time  integration  for  a  general  governing  equation  in  Eq.  (15)  follows  below. 

au  |  dF(u)  |  SFjvyv) 

dt  dx  dx  (15) 

Here,  U  is  solution  vector,  F  (U)  is  inviscid  flux  vector,  Fv  (U,  VU)  is  viscous  flux  vector,  and  G  (U)  is  source 
vector.  Approximation  of  solution  vector  U  is  written  in  terms  of  basis  function,  (j)1,  and  degree  of  freedoms  U]  as 
written  in  Eq.  (16).  Substituting  Eq.  (16)  in  Eq.  (15)  and  multiplying  resultant  equation  by  a  test  function  4>'j  gives  us 
a  weak  form.  Upon  integration  by  parts  to  both  the  inviscid  and  viscous  flux  terms,  the  weak  form  is  given  by  Eq. 
(17).  The  coefficient  of  the  first  term  in  Eq.  (17),  is  the  mass  matrix,  M  given  in  Eq.  (18). 

u^uj=±^]u,j 

I=°  (16) 
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jMM’ 


dx 


8{U‘ 


dt 


dx+  |  <f).F  |*  - 


dx+  |  <f>.Fv  |*=  J  dx 


(17) 


[M] 


At 


d(j) 

dx 


F\dx- 1  <j).F  |!  + N  —  .F  Life- 1  <t>.Fv  |!  +J  {(j).G)  dx 


dx 


(18) 


We  have  used  Legendre  polynomials  for  the  basis  functions,  which  leads  to  a  diagonal  mass  matrix.  Overall 
Jacobian  matrix  obtained  is  structured  in  block-wise  fashion  and  thus  is  amenable  to  parallelization  for  faster 
computations.  This  task  shall  be  taken  up  in  the  future  work.  For  approximation  of  numerical  fluxes  we  have  used 
Local  Lax  Friedrichs  (LLF)  flux  for  inviscid  fluxes  and  BR-1  (Bassi  Rebay  -  1)  for  viscous  fluxes  in  equations  19 
and  20.  As  a  standard  practice,  we  use  auxiliary  equations  to  deal  with  viscous  terms  in  DG  method. 

Hllf  (a,b)  =  ^(F(a)  +  F(b)-a(b-a)) 


a  =  max 

min(«,Z))<5<max(fl,Z)) 


1 


dF 

dU 


« 


(19) 


Hv(a,b)  =  -(Fv(a)  +  Fv(b )) 


(20) 


Figure  3.  Residual  convergence  plot  for  a  test  case  designed  for 
thermal  ablation  problem  and  time-step  of  1.0e-06  sec.  Plot 
shows  near  quadratic  convergence  of  the  scheme. 


Since  problem  is  non-linearly  coupled,  we 
use  Newton’s  method  to  solve  the  implicit 
system  and  fast  convergence  of  the  method 
depends  on  the  accurate  evaluation  of  Jacobian. 

Our  initial  attempt  was  at  getting  the  analytical 
expressions  of  the  Jacobians.  We  have  already 
demonstrated  complexity  and  limitations  that 
abound  analytical  approach  to  evaluate  jacobians 
in  [24]  and  [28].  To  summarize,  approach  of 
evaluating  analytical  expressions  is  very 
cumbersome,  and  has  convergence  and  stability 

problems  for  thermal  ablation.  It’s  also  not  generalized  and  typically  requires  a  lot  of  work  for  a  new  given  problem. 
For  any  new  problem,  new  analytical  expressions  need  to  be  found,  and  such  expressions  become  large  and  difficult 
to  handle.  This  approach  is  also  more  prone  to  errors.  Also,  Jacobian  evaluation  at  the  boundary  is  crucial,  e.g.  if  a 
variable  at  the  boundary  is  fixed,  then  while  calculating  Jacobian,  it  should  be  treated  as  constant  when 
differentiating  w.  r.  t.  other  variables;  and  differentiation  w.  r.  t.  the  fixed  variable  should  not  be  considered.  Due  to 
all  the  above  reasons  we  took  to  numerical  approach  to  evaluate  these  jacobians,  which  proved  to  be  much  more 
accurate,  effective  and  convenient  to  use.  This  approach  is  also  not  dependent  on  the  given  system  of  equations  and 
is  thus  general  in  application.  We  don’t  need  to  change  this  part  of  code  for  any  system  of  equations.  For  numerical 
approach  of  evaluating  jacobians,  we  perturb  a  solution  variable  and  then  calculate  the  change  in  flux  term  that  it 
produces.  By  applying  simple  finite  difference  we  get  a  first  order  approximation  to  the  Jacobian,  which  is  accurate 
enough  if  the  perturbation  is  sufficiently  small.  Due  to  coupling  between  elements  (through  derivatives  of  F+  and  F"), 
and  numerical  evaluation,  we  were  able  to  achieve  near  quadratic  convergence  for  thermal  ablation  problems  as 
shown  in  Fig.  3. 


B.  2-D  Discontinuous  Galerkin  method 
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For  the  time  integration  of  2-D  Discontinuous  Galerkin  method,  we  use  implicit  scheme.  2-D  model  problem  is 
given  by  Eq.  21,  where  F  and  Fv  have  x  and  y  components.  Multiplying  Eq.  21  by  a  test  function  and  integrating 
over  domain  (along  with  integration  by  parts  for  viscous  and  inviscid  flux  terms)  gives  us  weak  form  in  Eq.  22. 


dU 

dt 


+  V.K(U)-V.Fv(U,VU) 


=  G 


(21) 


d Cl  +  ^  vF.nd a  -  V vFdD.  -  [  ^  vFv .nd a  +  V v.Fvd Q  - 


vGd  Q  =  0 


(22) 


Basis  functions  chosen  for  the  above  formulation  are  Jacobi  polynomials,  and  are  given  in  Table  2  below.  These 
basis  functions  were  taken  from  [25].  Here,  coordinates  x  and  y  are  given  for  a  reference  square  element,  which  can 
be  transformed  to  a  general  quadrilateral  in  the  physical  space  using  bilinear  mapping.  These  basis  functions  are 
orthogonal,  and  thus  the  mass  matrix  will  be  diagonal,  but  since  we  allow  general  quadrilaterals  whose  opposite 
sides  may  not  be  parallel,  therefore  the  Jacobian  of  transformation  will  no  longer  be  independent  of  position  within 
the  reference  element,  and  hence  we  will  have  a  mass  matrix  that  will  no  longer  be  diagonal,  but  symmetric. 


Table  2:  Basis  functions  expressions  based  on  their  orders.  Here  x  and  y  are  coordinates  for  the 
reference  geometry  which  maps  to  a  generalized  quadrilateral. _ _ 


Basis/Order 

1 

2 

3 

4 

1 

1 

1 

1 

1 

2 

-l+2x 

-l+2x 

-l+2x 

3 

-l+2y 

-l+2y 

-l+2y 

4 

l-6x+6x‘ 

l-6x+6x3 

5 

l-2y-2x+2xy 

l-2y-2x+4xy 

6 

l-6y+6y~ 

l-6y+6y~ 

7 

-1+12x-30x2+20x3 

8 

- 1  +2  y+6x- 1 2  xy-  6y~+ 1 2  xy~ 

9 

- 1  +6  y-  6y~+2  x- 1 2  xy+ 1 2  yx~ 

10 

-l+12y-20y-+20y3 

IV.  Results  and  Discussions 

We  present  results  for  1-D  thermal  ablation  problem  that  was  described  in  section  111  and  also  demonstrate  our 
development  of  2-D  discontinuous  Galerkin  methods  and  its  application  to  a  double  mach  reflection  problem. 
Before  presenting  these  results  and  the  discussion  that  follows,  we  validate  our  DG  code  for  a  1-D  Sod’s  shock  tube 
problem,  comparing  results  with  the  analytical  solution. 

A.  Validation  case  for  discontinuous  Galerkin  method 

We  validate  results  for  a  shock  tube  problem  with  a  test  case  from  Toro  [29].  Sod’s  shock  tube  is  a  well  known 
problem  and  its  analytical  solution  is  also  available  for  test  of  accuracy  and  convergence.  Total  domain  length  is  1 
m,  and  is  divided  into  two  chambers,  left  and  right  at  x  =  0.5  m.  Initial  values  for  density  and  pressure  in  left 
chamber  are  both  1  in  SI  units,  and  in  right  chamber  are  0.125  and  0.1  respectively.  Velocity  is  initialized  to  a  zero 
value.  The  solutions  are  plotted  at  total  time  of  0.2  sec. 

This  case  is  run  using  implicit  DG  code  with  time-steps  of  1.0X10'3  and  l.OxlO"4  sec.  We  use  local  Lax- 
Friedrichs  for  numerical  fluxes,  and  show  comparison  of  solutions  of  order  of  accuracy  up  to  4  (LLF  01-04)  for  the 
code  with  the  exact  solution  in  Fig.  4.  Solution  accuracy  is  definitely  improved  at  a  lower  time-step.  Implicit  time- 
integration  scheme  has  also  has  its  own  dissipation  mechanism,  which  can  be  seen  in  smearing  of  solution  and 
oscillation  peaks  with  time  step  of  1.0X10'3  sec.  We  also  provide  close-ups  of  three  regions  namely,  expansion  fan. 
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a)  At  =  l.OxlO  03  sec  b)  At  =  l.OxlO04  sec 

Figure  4.  Validation  cases  for  shock  tube  problem  for  two  time  steps,  a)  1.0e-03  sec  and  b)  1.0e-04  sec. 
Solution  from  DG  simulations  for  high  order  accuracy  is  compared  to  exact  solution,  and  higher  accuracy 
gives  closer  solution. 


contact  discontinuity  and  shocks  to  give  a  better  idea  of  how  solution  approximations  improve  with  order  of 
accuracy  at  each  of  these  zones.  These  results  are  shown  in  Fig.  5  for  time-step  of  l.OxlO"4  sec. 


0.3 


0.25 


0.15 


0.1 


- LLF01 

-  LLF02 

-  LLF  03 

-  LLF  04 


a)  Expansion  wave  b)  Contact  discontinuity  c)  Shock  wave 

Figure  5.  Close  ups  of  three  regions  in  shock-tube,  a)  expansion  wave,  b)  contact  discontinuity  and  c) 
shock  wave.  Time  step  for  above  simulations  is  l.OxlO"4  sec.  Higher  polynomial  approximations  lead  to 
more  accurate  results. 


B.  1-D  thermal  ablation  using  discontinuous  Galerkin  method 

An  Arc -jet  wind  tunnel  experiment  was  conducted  in  [20]  with  a  carbon-phenolic  sample  that  was  exposed  to  a 
heat  flux  of  1400  W/cm2.  Experiment  results  were  compared  with  simulation  results  from  Charring  Material 
Ablation  code  or  CMA.  Later  on,  these  results  were  also  published  in  [1]  and  [2].  Total  run  time  for  this  problem  is 
5  sec,  and  surface  pressure  and  enthalpy  at  the  exit  are  0.22  atm  and  23300  J/g.  Total  thickness  of  the  test  sample  is 
1  cm. 

Here  we  consider  two  cases,  as  taken  in  [1]  and  [2].  Differences  in  both  these  cases  were  noted  in  the 
introduction  section.  To  compare  the  effect  of  diffusion  term  modeling  between  both  the  cases,  we  use  friction  term 
and  thermal  conductivity  from  [2]  in  both  the  cases.  First  we  present  results  obtained  for  the  case  with  D  =  0,  which 
was  considered  in  [1].  In  Fig.  6  a)  and  b),  we  show  spatial  distribution  of  solution  variables  at  a  time  of  4  sec.  It  can 
be  seen  from  plot  of  resin  density,  p,  that  pyrolysis  zone  extends  from  2  mm  to  4  mm.  Other  2  zones,  i.e.  virgin 
plastic  (pr  =  250  kg/m3)  and  char  zone  (p,  =  0  kg/m3)  can  also  be  identified  from  this  plot.  Pressure  varies  from  a 
value  of  nearly  8  atm  in  the  region  of  virgin  plastic,  where  it  is  flat,  to  a  low  0.22  atm  at  the  right  side  of  domain. 
High  pressure  gradient  of  the  gas  is  observed  within  the  pyrolysis  zone.  The  gas  flow  is  actually  driven  by  the 
gradient  of  the  product  of  the  gas  pressure,  P  and  void  fraction,  s.  This  is  seen  in  the  jump  of  velocity  magnitude  in 
the  pyrolysis  zone,  and  higher  velocities  are  attained  at  the  right  boundary  due  to  a  steeper  gradient  of  sP.  We  also 
observe  some  numerical  artifacts  in  the  solution  of  gas  momentum,  possibly  due  to  dominant  friction  term  in  Eq.  (3). 
Cooling  effect  of  pyrolysis  gas  can  also  be  seen  from  these  plots. 
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a)  Gas  density,  gas  momentum  and  pressure  b)  Gas  velocity,  resin  density  and  temperature 


c)  Temperature  distribution  at  different  times 

Figure  6.  Top  two  plots  a)  and  b)  show  all  solution  variables  for  Arc-jet  case  with  D  =  0  at  time  =  4  sec.  a) 
shows  gas  density,  spg,  gas  momentum,  epgu  and  pressure,  P;  b)  shows  gas  velocity,  u,  resin  density,  pr,  and 
temperature,  T;  and  c)  shows  temperature  distribution  within  the  domain  at  times  =  0.1, 1,  2,  3,  4,  and  5  sec. 

Temperature  distribution  within  the  material  has  steeper  gradient  at  the  pyrolysis  front,  seen  at  around  2  mm  in 
Fig.  6b),  as  compared  to  elsewhere  in  domain.  There  are  two  heating  mechanisms  in  the  domain,  thermal  conduction 
and  convection  cooling.  In  beginning  of  pyrolysis  zone,  we  expect  thermal  conduction  to  be  the  dominant 
phenomena,  whereas  at  other  places  in  the  char  we  expect  both  convection  cooling  and  conduction  to  play  an 
equally  important  role.  This  explains  the  observed  temperature  distribution  and  the  trend  of  thermal  response  (at 
different  times)  that  can  be  seen  in  Fig.  6c),  where  a  steep  rise  in  temperature  is  slowed  down  by  the  convection 
cooling  of  the  exiting  pyrolysis  gas.  Based  on  our  calculations,  we  find  that  temporal  increase  in  pyrolysis  gas 
density  in  the  virgin  plastic  zone  (region  to  the  left  of  pyrolysis  zone)  is  due  to  resin  decomposition  rate,  R,  which 
also  leads  to  corresponding  gas  generation.  The  temperature  in  the  virgin  zone  is  300  K,  so  value  of  R  should  be 
very  close  to  zero,  which  is  not  the  case  with  the  model  for  R  given  in  [1]  and  [2].  Based  on  these  calculations,  we 
come  to  conclusion  that  all  the  gas  that  is  generated  during  pyrolysis  (in  the  pyrolysis  zone)  primarily  leaves  out  to 
the  atmosphere.  Both  friction  and  inertia  terms  in  [2]  provide  resistance  to  the  flow  of  gas.  Fligher  values  for  friction 
term  in  [1]  results  in  higher  resistance  to  the  flow  and  hence  larger  gas  density  values  are  observed  within  the  virgin 
plastic  domain.  This  would  result  in  higher  pressure  to  be  observed  for  [1]  as  compared  to  [2].  This  difference  was 
noted  in  the  pressure  results  reported  in  [1]  and  [2].  Gas  pressure  for  the  day  probe  in  the  virgin  plastic  domain  is 
nearly  75.67  atm  for  [1]  and  30  atm  for  [2]  at  a  flight  time  of  27  sec.  We  are  investigating  this  difference  in  two 
models  between  2  papers.  We  will  note  here  that  we  have  used  friction  model  from  [2]  for  both  the  cases  that  we 
ran,  to  compare  the  effect  of  Diffusion  term  alone  between  both  the  cases. 

In  Fig.  7,  we  compare  our  results  for  temporal  variation  of  surface  temperature  (on  the  right  side  of  domain 
shown  in  Fig.  1).  We  provide  results  for  our  simulation  for  both  curve  fit  -  1  and  curve  fit  -  2  for  internal  energy,  eg. 
As  we  have  already  mentioned  that  curve  fit  -  2  is  more  accurate  than  curve  1,  we  find  that  our  results  match  the 
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Figure  7.  Temporal  variation  of  surface  temperature  for  DG  simulation  with  D  =  0,  and  comparison  with 
1998  paper  and  experimental  results. 

experimental  results  for  initial  1  sec,  and  after  that  due  to  high  gas  velocity,  there’s  no  further  rise  in  temperature 
observed.  High  gas  velocity  leads  to  high  cooling  at  the  surface,  which  exactly  balances  out  net  heat  influx  to  the 
surface  and  thermal  conduction. 

Second  case  that  we  present  is  with  nonzero  diffusion  term  (D^0)  in  Eq.  2.  This  case  is  exactly  similar  to  [2]. 
Shown  below  in  Fig.  8,  are  results  for  solution  variables  at  time  of  4  seconds.  We  now  highlight  the  differences  in 
both  the  cases,  with  and  without  diffusion  term.  Diffusion  term  in  momentum  equation  is  related  to  the  second 
derivative  of  gas  pressure,  and  results  in  more  diffused  profile  for  gas  density.  Pressure  is  related  to  gas  density  and 
temperature,  and  both  are  more  diffused  out,  which  causes  a  lower  pressure  (actually  eP)  gradient  that  drives  the 
flow,  and  therefore  we  get  lower  gas  momentum  near  exit,  as  compared  to  case  D  =  0.  Therefore  we  find  these 
results  are  consistent  with  basic  gas  dynamics. 

Due  to  diffusion  term,  gas  velocity  at  the  point  of  gas  generation  is  lower,  and  hence  gas  density  is  higher 
compared  to  case  with  D  =  0  in  the  pyrolysis  zone.  Larger  gas  density  absorbs  more  heat  from  the  material,  and 
hence  temperatures  at  pyrolysis  zone  are  lower  for  this  case.  This  in  turn  affects  the  decomposition  of  resin  material 
and  we  see  a  broader  pyrolysis  zone  in  this  case,  which  extends  roughly  from  3  mm  to  7  mm.  This  diffusion 
mechanism  is  in  general  better  for  the  vehicle’s  body  since  it  provides  for  more  time  that  the  vehicle  gets  before  its 
temperature  starts  to  increase.  For  example,  temperatures  above  500  K  are  observed  at  location  of  2.25  mm  for 
earlier  case  (D  =  0)  and  for  this  case  at  location  of  3.17  mm  at  time  of  4  sec.  There’s  no  significant  difference  in 
pressure  for  both  cases,  pressure  in  this  case  is  approximately  0.5  atm  higher  than  that  of  case  -  1.  Velocities  in  the 
domain  for  this  case  are  lower  as  compared  to  case  with  D  =  0,  therefore  we  have  lower  cooling  within  the  domain, 
and  since  heat  conduction  is  a  diffusion  mechanism,  we  observe  more  diffused  profiles  of  temperature.  So,  from 
looking  at  temperature  profile  we  can  tell  which  effect  is  dominant,  conduction  or  convective  cooling. 

We  also  show  comparison  of  results  obtained  for  both  D  =  0  and  D  ^  0  at  final  time  of  5  sec  below  in  Fig.  10, 
and  we  can  easily  see  that  lesser  amount  of  resin  material  is  used  for  the  case  with  D  ^  0  (top  left  plot  in  Fig.  10). 
The  pyrolysis  zone  is  identified  from  0.002  m  to  0.006  m  for  non-zero  diffusion  case,  and  that  for  zero  diffusion 
from  0.0005  m  to  0.0025  m;  therefore  pyrolysis  zone  is  twice  wide  for  non-zero  diffusion  case.  Also,  the  pyrolysis 
front  for  non-zero  diffusion  case  travels  at  lower  speed  than  case  with  zero  diffusion.  Since  temperatures  around  the 
pyrolysis  zone  are  lower  for  non-zero  diffusion  than  temperatures  around  pyrolysis  zone  for  zero  diffusion,  net 
decomposition  rate  for  non-zero  diffusion  case  is  lower,  hence  speed  of  pyrolysis  front  is  slower.  Lower  speeds  for 
exiting  gas  within  the  domain  can  be  seen  in  Fig.  10  for  non-zero  diffusion  and  hence  larger  gas  densities,  which 
means  more  absorption  of  heat  from  the  material  and  that  is  the  reason  why  we  observe  lower  and  more  diffused 
temperature  plots. 

Temporal  variation  of  surface  temperature  with  simulation  results  by  [2]  and  experimental  results  by  [20]  are 
given  in  Fig.  9.  Again  we  plot  the  results  for  both  curve  fits  -  1  &  2  for  internal  energy.  Our  surface  temperature 
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shoots  to  higher  value  at  later  time  than  experimental  results.  Reason  for  this  as  already  explained  is  lower  cooling 
for  this  case.  (Given  our  lack  of  clarity  about  some  aspects  of  their  model,  our  plot  is  close  to  experiment  results) 


a)  Gas  density,  gas  momentum  and  pressure  b)  Gas  velocity,  resin  density  and  temperature 


c)  Temperature  distribution  at  different  times 

Figure  8.  Top  two  plots  a)  and  b)  show  all  solution  variables  for  Arc-jet  case  with  non-zero  diffusion  term  at  time 
=  4  sec.  a)  shows  gas  density,  spg,  gas  momentum,  ep„u  and  pressure,  P;  b)  shows  gas  velocity,  u,  resin  density,  pr, 
and  temperature,  T;  and  c)  shows  temperature  distribution  within  the  domain  at  times  =  0.1, 1,  2,  3,  4,  and  5  sec. 
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Figure  9.  Temporal  surface  temperature  variation  for 
non  zero  D  and  2002  paper 


c)  Gas  density  d)  Gas  velocity 

Figure  10.  Comparison  plots  for  cases  with  zero  and  non-zero  diffusion  coefficients  for  a)  resin  density,  b) 
temperature,  c)  gas  density  and  d)  gas  velocity  at  time,  t  =  5  sec.  We  can  notice  for  same  heat  in-flux  and  total 
time,  lesser  material  is  used  with  non-zero  diffusion  case  than  with  zero  diffusion  case. 
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Figure  11.  Schematic  of 
problem  (taken  from  [26]) 


2-D  double  mach  reflection 


C.  2-D  Discontinuous  Galerkin  application  to  Double  Mach  Reflection  Problem 

In  our  attempt  to  develop  2-D  capability  for  discontinuity  Galerkin  method,  we  applied  this  code  to  problems  of 
advection  and  Burgers  equations  and  double  mach  reflection  problem.  For  time-integration  we  have  used  implicit 
method  (BDF  -1  or  Backward  Difference  Formula). 

Here  we  only  present  results  for  double  mach 
reflection.  This  is  a  well  known  problem,  and  has 
been  attempted  analytically,  experimentally  and 
numerically.  In  the  field  of  DG  methods,  it  was 
solved  in  [6].  We  try  this  problem,  as  a  first  step  into 
building  a  2-D  DG  capability  in  our  in-house  code, 

MIG.  Schematic  of  this  problem  (taken  from  [26]) 
is  given  in  Fig.  11. 

Here  supersonic  flow  meets  the  ramp  (horizontal 
line  in  Fig.  11)  at  an  angle  of  a,  due  to  which  the 
incoming  shock  is  reflected  from  this  compressive 
ramp,  and  this  interaction  is  referred  to  as  double 
mach  reflection.  For  boundary  conditions,  and 
domain  information  please  refer  to  [26].  Since  this 
problem  has  a  strong  shock  and  due  to  absence  of 
an  external  artificial  dissipation  mechanism  in  our 
scheme,  we  get  overshoots  and  undershoots  across 
the  shock;  this  causes  the  density  and  total  energy  to 
take  negative  values  near  shock.  In  our  code,  we 
limit  values  of  density  and  total  energy  to  their 
minimum  physically  possible  values.  By  this  we 
have  been  able  to  get  results  for  k  =  1  (linear 
polynomial  approximation)  for  DG  method.  For 
higher  k  (>=  2),  the  jumps  across  shock  become 
steeper  and  we  get  convergence  issues,  so  we 
haven’t  been  able  to  get  solutions  for  full  run  time 
for  higher  k.  In  order  to  resolve  these  issues,  we  are 
working  on  implementing  artificial  dissipation 
scheme  for  removing  these  oscillations  and  as 
shown  in  [12],  with  higher  k,  we  should  expect 
more  accurate  solutions. 

Fig.  12  shows  density  contours  for  gas  density  at 
t  =  0.2  sec.  For  these  simulations  we  use  time-step 
of  1.0e-4  sec,  and  the  results  are  shown  for  40  X  10, 

80  X  20,  and  160  X  40  mesh  and  k  =  1,  i.e.  linear 
polynomial  approximation.  The  domain 

dimensions  are  4  m  in  x-direction  and  1  m  in  y- 
direction.  Initial  location  of  shock  is  at  (1/6,  0),  and 
aligned  at  60°  to  x-axis.  We  can  easily  notice  that 
solution  quality  clearly  improves  with  refined  mesh, 
but  due  to  lack  of  artificial  dissipation  for  this 
problem,  the  jumps  across  the  shock  become  more 
severe  and  hence  solution  becomes  difficult  to 

converge.  The  solution  contours  are  comparable  to  other  published  results  of  this  problem,  but  with  higher  mesh  we 
can  observe  oscillatory  nature  of  density  contours.  At  very  coarse  mesh,  this  is  not  observed,  but  many  fine  flow 
features  become  visible  with  finer  mesh.  Also  with  k  =  2,  i.e.  quadratic  polynomial  approximation  jumps  become 
more  severe  and  simulations  hard  to  converge. 

Convergence  results  for  double  mach  reflection  problem  are  also  shown  in  figure  given  below,  in  Fig.  13.  Data 
points  for  both  L2  and  Loo  norm  of  density  are  plotted  for  all  the  three  meshes  considered  above,  along  with  their 
linear  fit,  and  it  can  be  seen  easily  that  the  simulations  converge  to  a  fixed  solution  with  the  given  three  meshes. 
Norms  are  plotted  against  mesh  measure,  which  is  the  uniform  edge  thickness  of  the  meshes  used,  which  are  same 
for  x  and  y  axis. 


Figure  12.  Density  contours  obtained  for  40  X  10,  80  X 
20,  and  160  X  40  mesh  (arranged  from  top  to  bottom) 
and  k  =  1  (linear  polynomial  approximations),  time-step 
for  above  simulations  is  1.0e-04  sec,  and  results  are 
shown  at  t  =  0.2  sec 
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Figure  13.  Mesh  refinement  study  for  the  double  mach  reflection  problem, 
where  norm  of  the  density  is  plotted  against  the  mesh  measure,  h 


V.  Conclusion  &  Future  Work 

An  implicit  discontinuous  Galerkin  method  to  solve  thermal  ablation  is  introduced  and  the  results  are 
benchmarked  with  the  previously  reported  results.  We  are  able  to  achieve  near  quadratic  convergence  using 
numerically  evaluated  jacobians  in  Newton  Raphson  scheme.  Two  cases  with  and  without  the  diffusion  term  were 
compared  and  it  was  found  that  the  modeling  of  the  dynamics  of  pyrolysis  gas  does  affect  the  net  thermal  response 
of  the  TPS.  The  case  with  D  ^  0  has  more  diffused  profiles  for  pressure,  density  and  temperature  and  hence  the 
pyrolysis  front  travels  at  a  slower  speed  in  this  case  than  for  D  =  0.  Consequently,  lower  amount  of  resin  material  is 
used  up  with  D  ^  0.  Hence,  neglecting  the  diffusion  term  will  lead  to  over  prediction  of  TPS  mass,  which  means 
increased  fuel  cost  for  the  space  mission.  So,  we  see  that  modeling  of  motion  of  pyrolysis  gas  is  very  crucial  in 
accurately  deciding  the  net  response  of  TPS,  and  this  is  neglected  in  many  numerical  codes  on  in-depth  thermal 
ablation.  Pyrolysis  gas  was  considered  to  be  in  thermal  equilibrium  in  current  work,  and  we  will  explore  avenues  for 
considering  finite  rate  reaction  in  our  future  work.  Some  modeling  inconsistencies  in  between  [1]  and  [2]  are  also 
noted  and  this  will  be  further  investigated.  To  extend  our  current  work  to  higher  dimensions,  we  are  building  a  2-D 
DG  capability,  and  as  a  part  of  our  ongoing  work  in  this  area,  we  have  shown  results  for  double  mach  reflection 
problem.  As  a  next  step,  we  plan  to  apply  artificial  dissipation  to  this  2-D  double  mach  reflection  problem,  which 
will  solve  the  problem  of  capturing  shock,  without  spurious  oscillations  of  jumps  (overshoots  and  undershoots  at 
discontinuities). 
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